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Abstract 

We present results from a detailed Quantum Monte Carlo 
study of BEC applied to JILA experiment [Jin et al, Phys. Rev. 
Lett.78,764,1997][l]. This is the first Monte Carlo approach ( based 
on Feynman-Kac path integral method) to the above problem where 
good qualitative agreement is found for both the lowest lying m = 2 
and m = mode. We found an upward shift of the experimental data 
for m = mode at around T = 0.7Tq (Tq is defined as the predicted 
BEC transition temperature for a harmonically confined ideal gas) 
when the effect of noncondensate was considered. 
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1 Introduction 



After the experimental realization of Bose Einstein Condensation in alkali 
vapors in 1995 [2], and subsequent experiments pertinent to temperature de- 
pendence of frequencies and damping rate[l], there have been a lot of the- 
oretical studies [3-9] to explain the experimental observations in connection 
with temperature dependent frequency shifts corresponding to different an- 
gular momenta, m=0 and m=2 modes in particular [2]. Results have been 
reported in which theoretical data agreed well with experimental values for 
m=0 mode showing an upward trend of frequencies with rise in temperature 
[4,9]. But in all cases the agreement is rather poor when it comes to m=2 
mode. There is an agreement up to T = 0.6T , beyond which frequencies rise 
with increase in temperature deviating from the downward trend of experi- 
mental data. In this article, we would like to report a diffusion Monte Carlo 
study of the frequency shifts of m=2 and m=0 modes in a dilute gas of Rb 87 . 
In our non mean field study, we see agreement with experimental study (Fig 
3a of Ref 2]for m = 2 mode all the way to T = 0.9T (Fig 6). When we 
consider the the dynamics of the thermal cloud separately, the upward shift 
( Fig 7) at T = 0.7T which is similar to JILA experiment is observed for 
m = mode. This agrees with the results obtained from the revised gapless 
theory of Morgan [6, 7]. 

The dynamical behavior of dilute alkali BECs at T=0 can be well 
described by Gross-Pitaevskii eqn(GPE)[10]. 

in W££ = [-^-| + vUf) + gW,t)\ 2 W,t) (l) 

where g = 4**-2, V is the scattering length and 'm^' is the mass of Rb atom. 
But it seems to be inadequete at finite temperatures. The total density of 
the atoms is related to BEC density and normal component as follows: At 
T=0 the normal component is not equal to zero in the interacting case and 
is referred to as 'quantum depletion'. At finite temperature, thermal atoms 
also contribute to the normal component. Since mean field wavefunctions do 
not account for the normal component, it gives accurate energy spectrum if 
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depletion is smallfll]. Near To, the quantum depletion becomes significant 
and mean field treatment breaks down. 

The other mean field theories which have been used so far, are 
based on HF and HFB-Popov equations [12] and break down near T as 
its effective single particle spectrum always displays a gap. In 1998, a self 
consistent gapless non-divergent theory [5] was developed and a closed set of 
coupled equations were solved numerically. In this analysis, only the dynam- 
ics of condensate was considered and downshift of data was observed for both 
m = and m = 2 mode. Subsequently, with the more sophisticated theory 
of Morgan[6] an upward shift at T = 0.6T was achieved[7,8] Analytic ex- 
pressions [13] for temperature dependent frequencies were obtained for m=0 
and m=2 by time dependent variational technique. 

Eventhough Ref[7] has the best agreement with JILA data till 
date, solving coupled partial differential equations numerically is not an easy 
task. The chief purpose of this paper is to go beyond mean field theory 
with a comparatively simpler numerical procedure which would work at all 
temperatures. We propose to explore finite temperature aspect of BEC by 
quantum nonperturbative technique, namely Feynman- Kac (FK) [14-16] pro- 
cedure. To be precise, we use Generalised Feynman-Kac ( GFK ) method[17] 
to make the rate of convergence faster. Since Quantum Monte Carlo meth- 
ods are computationally expensive, we are simulating only 2000 interacting 
atoms at this moment. Increasing number of interacting atoms would change 
our results quantitativelynot qualitatively. This paper is organized as fol- 
lows : In Sec 2, we discuss the path integral technique at zero and finite 
temperature as a many body techinique, the Schroedinger formulation of Rb 
condensate and noncondensate, fundamental concepts of BEC and finite tem- 
perature excitations. In Sec 3, we discuss the numerical procedure. In Sec 4, 
we present all the numerical results pertinent to energies and frequencies at 
different temperature. Finally in Sec 5, we summarize our results. 
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2 Theory 



To connect Feynman-Kac or Generalized Feynman Kac ( GFK ) to other 
many body techniques our numerical procedure ( GFK ) [18-19] has a straight- 
forward implementation to Schroedinger's wave mechanics. Since at low tem- 
perature the de Broglie wavelength of the atoms become appreciable, we do 
a full quantum treatment. GFK is essentially a path integral technique with 
trial functions for which operations of the group of the wave function keep 
points in the chosen nodal region, provide an upper bound for the lowest state 
energy of that symmetry. The nodal region with the lowest energy serves as a 
least upper bound. If the nodal region has exact nodal structures of the true 
wave function the random walk is exact in the limit scale, time for walk, and 
number of walks get arbitrarily large. To calculate energy we approximate 
an exact solution ( i.e., the GFK representation of it ) to the Schroedinger's 
equation, whereas most of the other numerical procedures approximate a so- 
lution to an approximate Schroedinger equation. From the equivalence of 
the imaginary time propagator and temperature dependent density matrix, 
finite temperature results can be obtained from the same zero temperature 
code by running it for finite time. So from all these aspects, Generalised 
Feynman-Kac method turns out to be a potentially good candidate as a 
sampling procedure for Bose gases at all temperatures. Next we consider the 
Feynman-Kac formalism and then show how it can be modified to get the 
Generalized Feynman-Kac version of it. 
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2.1 Path integral Theory at T=0 
2.1.1 Feynman-Kac Path integretion 

For the Hamiltonian H = —A/2 + V(x) consider the initial value problem 

= (4 + i> (t ,*) 

u(0,x) = f(x) (2) 

with x G R d and u(0, x) = 1. The solution of the above equation can be 
written in Feynman-Kac representation as 

u(t, x) = E x exp[- f V(X(s))ds] (3) 
Jo 

where X(t) is a Brownian motion trajectory and E x is the average value of 
the exponential term with respect to these trajectories. The lowest energy 
eigenvalue for a given symmetry can be obtained from the large deviation 
priniciple of Donsker and Varadhan [20] , 

li = - lim -ln\E x exp\- f V(X(s))ds\] (4) 

t-HX> t JO 

The above formalism is valid for any arbitrary dimension d ( for a system 
of N particles in three dimensions d = 3N). Generalizations of the class of 
potential functions for which Eqns. 3 and 4 are valid are given by Simon[21] 
and include most physically interesting potentials, positive or negative, in- 
cluding, in particular, potentials with 1/x singularities. It can be argued that 
the functions determined by Eq(3) will be the one with lowest energy of all 
possible functions independent of symmetry. Restrictions on allowed Brow- 
nian motions must be imposed to get a solution of the desired symmetry if 
it is not the lowest energy solution for a given Hamiltonian. Since the above 
energy formula gives the lowest energy corresponding to any symmetry, the 
same formula can be used to calculate ground and excited states of a quantum 
mechanical system. Although other interpretations are interesting, the sim- 
plest is that the Brownian motion distribution is just a useful mathematical 
construction which allows one to extract the physically relevant quantities, 
the ground and excited state energy of a quantum mechanical system. In 
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numerical implementation of Eq(4) the 3N dimensional Brownian motion is 
replaced by 3N independent, properly scaled one dimensional random walks 
as follows. For a given time t and integers n and 1 define [18] the vector in 

W(l) = W(t,n,l) = (w 1 1 (t,n,l),w 2 1 (t,n,l),w 3 1 (t,n,l).... (5) 

Wi N (t, n, l)w 2 N (t, n, l)w 3 N (t, n, I) 

where 

<(W) = E% ( 6 ) 

k=l v n 

with w/(0, n,l)=0 for i = 1, 2, N;j = 1, 2, 3 and 1 = 1,2, ,nt. Here 

e is chosen independently and randomly with probability P for all i,j,k such 
that P(e % j k = l)=P(e' t jk = — 1)—\- It is known by an invariance principle[22] 
that for every v and W(l) defined in Eq(5) 

-i nt 

n ¥/(;EW0))^ (7) 
u i=i 
t 

= P(J V(X(s))ds < v 
o 

Consequently for large n, 

t 

P[exp(- J V{X{s))ds) < is] (8) 
o 

i nt 

w P[exp( Y,V(W(l))) < v\ 

n i=i 

By generating N rep independent replications Zi,Z 2 ,....Z Nrep of 



i nt 

Z m = eM-(—Y / V(W(l))) (9) 
n i=i 

and using the law of large numbers, [Z\ + Z 2 + ...ZN rep )/N rep = Z{t) is an 
approximation to Eq(3) 

H » ~logZ(t) (10) 

Here W m (l),m = l,2,N rep denotes the m th realization of W(l) out of N rep 
independently run simulations. In the limit of large t and N rep this ap- 
proximation approaches an equality, and forms the basis of a computational 
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scheme for the lowest energy of a many particle system with a prescribed 
symmetry. In dimensions higher than 2, the trajectory x(t) escapes to infin- 
ity with probability 1 . As a result, the important regions of the potential 
are sampled less and less frequently and the above equation converges slowly. 
Now to speed up the convergence we use Generalized Feynman-Kac (GFK) 
method. 

2.1.2 Generalized Feynman Kac path integretion 

To formulate the generalized Feynman-Kac method we first rewrite the 
Hamiltonian as H = H + V p , where H = —A/2 + /j, T + At/j T /2ip T and 
V p = V — {jiT + AiPt/Z^t)- Here ipx is a twice differentiate nonnegative 
reference function and Hipx = \it^t- The expression for the energy can now 
be written as 

fi = fi T ~ hm -ln[E x exp[- f V p (Y(s))ds}} (11) 

t^oo t JO 

where Y(t) is the diffusion process which solves the stochastic differential 
equation 

mt)= ^§w dt+mt) (12) 

The presence of both drift and diffusion terms in this expression enables the 
trajectory Y(t) to be highly localized. As a result, the important regions of 
the potential are frequently sampled and Eq (11) converges rapidly. 



2.2 Path integral theory at finite temperature 

The temperature dependence comes from the realization that the imaginary 

time propagator k(2, 1) is identical to the temperature dependent density 

matrix p(2, 1) if t =>- /3 = 1/T holds. 

This becomes obvious when we consider the eqs[23] 

dk(2, 1) 
dt 2 

and 



H 2 k(2,l) (13) 



■0 = ^(2,1) (14) 



and compare 

fc(2,l) =Y.H^W{x l )e-^- t ^ (15) 

i 

and 

p(2,l) = ^0 i (x 2 )^(x 1 )e-^ (16) 

i 

For Zero temp FK we had to extrapolate to t =>- oo. For finite run time t in 
the simulation, we have finite temparature results. In this section we show 
how we change our formalism to go from zero to finite temperature. We 
begin with the definition of finite temperature. A particular temperature 'T' 
is said to be finite if AE < kT holds. The temperature dependent density 
matrix can be written in the following form 

p(x,x',P) = p {0 \x,x',(3)x < exp[- j P V p [X{s)\ds\ > DRW (17) 

J o 

The partition function can be recovered from the above as follows: 

J p(x,x,/3)dx = J p {0 \x,x,/3)dxx < exp[- V p [X(s)\ds\ > D rw (18) 

In the usual notation, the above equation reads as 

Z(x,(3) = Z°(x,P)x < exp[- ( P V p [X{s)]ds\ > DRW (19) 

J o 

At finite temperature thus free energy can be written as 

F = -lnZ(x, (3) /I3 = -lnZ°(x, {3)/f3 - In < exp[- [ V p [X(s)]ds] >drw/P 

Jo 

(20) 
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2.3 Schroedinger Formalism for Rb condensate at 
T=0 

In the JILA experiment different frequency modes are labeled by their angu- 
lar momentum projection on the trap axis. In cylindrical symmetry, m = 2 
mode is an uncoupled one and there are two coupled oscillations for m = 
mode [24]. As a matter of fact Stringari [25] showed that m = mode is a 
coupled oscillation of a quadrupolar surface oscillation and a monopole. In 
the noninteracting case, these two modes are degenerate with uj/uj x = 2 . 
We choose to work in the cylindrical coordinates as the original experiment 
had an axial symmetry, cylindrical coordinates are the natural choices for 
this problem. We consider a cloud of N atoms interacting through repulsive 
potential placed in a three dimensional harmonic oscillator potential. At low 
energy the stationary state for the condensate can be represented as 

[-A/2 + V int + V trap ]%{f) = /i c ^ (f) (21) 



1 N 

[-A/2 + V mt + - ]>>, 2 + y 2 + \ a z 2 ]]M?) = ^o(r) (22) 
z i=i 

where \Y J f=i\ x i 2 + Di 2 + ^aZi 2 } is the anisotropic potential with anisotropy 
factor \ a = ^. Now 

V int = V Morse = ^Vinj) ^D^^^^ ~ 2 )] ( 23 ) 

In the above potential V ' is the location of the well minimum and 
V is the width of the Morse potential. The above Hamiltonian is not separa- 
ble in spherical polar coordinates beacause of the anisotropy. In cylindrical 
coordinates the noninteracting part behaves as a system of noninteracting 
harmonic oscillators and can be writtem as follows : 

\- — — ( 9 1 ° 2 1 92 

+ ^(p 2 + \a 2 Z 2 )}MP,z) 

= HcMPiz) (24) 
9 



The energy '//' of the above equation can be calculated exactly which is 

{2n p + \m\ + 1) + (n z + 1/2)A (25) 

In our guided random walk we use the noninteracting solution of Schroedinger 
equation as the trial function as follows [26] : 

^„, B , m (f) ^ exp^ H nz (z) x e^^Vf/' 2 ) (26) 

2.4 The effect of noncondensate 

In the case of noncondensate the system can be considered as a thermal gas. 
To calculate noncondensate energy and density we need to study the effect of 
noncondensate explicitly and consider the following stationary state for the 
thermal gas. 

[- A/2 + 2V int + Vtraptyj (f) = (f) (27) 

1 N 

[-A/2 + 2V int + - ]T [x 2 + y 2 + A^%](f) = Linear) (28) 
1 i=i 

The basis wavefunction ipj which describes the noncondensate should be 
chosen in such a way that it is orthogonal to tp as in Eq.(ll) The most 
common way to achieve a orthogonal basis in Schroedinger prescription is 
to consider the dynamics of noncondensate in an effective potential[6,27] 
V e ff = Vtrap + 2 Vi n t ■ The factor 2 represents the exchange term between two 
atoms in two different states. The energy in the case of lowest lying modes 
then corresponds to \x = fi c + fi nc . One can calculate the [i nc using the same 
parameters as discussed in Sec 3.1. 

2.5 Fundamentals of BEC 

Even though the phase of Rb vapors at T=0 is certainly solid, Bose con- 
densates are preferred in the gasous form over the liquids and solids because 
at those higher densities interactions are complicated and hard to deal with 
on an elementary level. They are kept metastable by maintaining a very 
low density. For alkali metals, 77, the ratio of zero point energy and molec- 
ular binding energy lies between 10~ 5 and 10" 3 . According to the theory of 
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corresponding states [28] since for the T=0 state of alkali metals, r\ exceeds 
a critical value 0.46, the molecular binding energy dominates over the zero 
point motion and they condense to solid phase. But again the life time of a 
gas is limited by three body recombination rate which is proportional to the 
sqaure of the atomic density. It gets suppressed at low density. Magnetically 
trapped alkali vapors can be metastable depending on their densities and life- 
times. So keeping the density low only two body collisions are allowed as a 
result of which dilute gas approximation [29] still holds for condensates which 
tantamounts to saying na 3 << 1 (a is the scattering length of s wave). Now 
defining n = N/V = r~ 3 as a mean distance between the atoms ( definition 
valid for any temperature ), the dilute gas condition reads as a « r av and 
zero point energy dominates (dilute limit). In the dense limit, for a ps r av 
on the other hand the interatomic potential dominates. The gas phase is 
accomplished by reducing the material density through evaporative cooling. 

2.6 Finite temperature Excitations : 

Finite temperature excitation spectrum is obtained by using the path inte- 
gral formalism used in Section 2.2. In our analysis, we first assume that the 
condensate oscillates in a static thermal cloud. There are no interactions 
between the condensate and the thermal cloud. The principal effect of finite 
temperature on the excitations is the depletion of condensate atoms. We 
want to calculate the collective excitations of Bose Einstein condensates cor- 
responding to JILA Top experiment ( m=2 and m=0 mode). Eventually for 
m = mode , we consider the effect of thermal cloud separately. 
Condensation fraction and Critical temperature : In the noninter- 
acting case for a harmonic type external force the theoretical prediction for 
condensation fraction is 




(29) 



Critical temperature can be defined as 



T = 



0.94 x TiuN 1 ' 3 
k B 



(30) 
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Q = (wfu;,) 1 ' 3 (31) 

From Eq.(29), we see that as temperature increases, condensation fracton 
decreases in the noninteracting case. Interaction lowers the condensation 
fraction for repulsive potentials. Some particles always leave the trap be- 
cause of the repulsive nature of the potential and moreover, if temperature 
is increased further, more particles will fall out of the trap and get thermally 
distributed. This decrease in condensation fraction eventually would cause 
the shifts in the critical temperature. We would observe this in Section 4.2 
(Fig. 4 ). Earlier this was done by W. Krauth[30] for a large number of 
atoms by path integral Monte Carlo method. In our analysis, we denote 'T ' 
as transition temperature following Ref[2] 
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3 Numerical procedure 



3.1 Dilute limit 

In the dilute limit and at very low energy only binary collisions are possible 
and no three body recombination is allowed. In such two body scattering at 
low energy first order Born approximation is applicable and the interaction 
strength 'D' can be related to the single tunable parameter of this problem, 
the s-wave scattering length 'a' through the relation given below. This single 
parameter can specify the interaction completely without the details of the 
potential in the case of pseudopotentials. We use Morse potential because 
it has a more realistic feature of having a repulsive core at r\j = than 
other model potentials. Secondly, using this realistic potential allows us to 
calculate the energy spectrum exactly as opposed to the case of 5 function 
potential where it is calculated perturbatively[31]. In our case the interaction 
strength depends on two more additional parameters, r and a. 

mD 



4nh 2 



V{r-)d A r (32) 



It is worth mentioning over here that instead of actual scattering length we 
use the Born approximation to it. Since we are dealing with a case of low 
energy and low temperature it is quite legitimate to use the above expression 
as a trickery to calculate the strength of Morse interaction [32]. As a matter 
of fact in Ref[33] the author has justified using Eq.(32) for a 5 function 
potential. So if it is justified to do it for 5 function potental it is even more 
justified to do so for Morse potential which is finite and short-ranged. 

The Morse potential for dimer of rubidium can be defined as 

J2V(rij) = Y J D[e~ ai - r ' ro \e' a{r - ro) - 2)] (33) 

where a is the depth of the Morse potential. Using the above potential 

D = 4 f W - (34) 

me aro (e aro - 16) 

The Hamiltonian for Rb gas with an asymmetric trapping potential and 
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Morse type mutual interaction can be written as 

N 

[-/> 2 /2m£vf + £y (r y 

i=l i,j 
^ JV N N 

+y K 2 E ^ 2 + < E tf 2 + ^ E *?) WO 

Z i=l i=l i=l 

= ^(r) (35) 

The above Hamiltonian can be rescaled by substituting r 1 = sr and \i = fx U 

as 

h 2 N Ah 2 nn^ 
\ V V- 2 + V [ e -a(»-5-ro)f e -a(ri5-n)) _ 2 )1 

L 2ms 2 ^ * ^m^e^fe" - 16) 1 V ;J 

ms 2 N N N 

+^K 2 E ^ 2 + ^ 2 E ^ 2 + ^ 2 E ^ 2 )]^ (r) 

^ 1=1 i=l 1=1 

= ^^(f) (36) 

[-V V, 2 -4 yr e -«(r«-ro)( e -a(r tf -ro) _ 2) 1 

2 ^ se ar °(e ar ° - 16) v ;J 

m 2 u; x 2 s 4: , . a;,, 2 . a;. ° 



-^^E(^ 2 + ^ 2 + ^ 2 M0 

= -^ *^-^(f) (37) 

Let = i =>. s 2 = _h_ ■ the natural unit of i en gth. Let = 1 =>• 

[/ = -^o = hw x is the natural unit of energy. Then the standard form of the 
equation becomes 

[ *k l h (e--16)^ ie (e 2)J 

-^E(^ 2 + ^ 2 + ^W) 

= -to^lf) (38) 
With uj x = oj y = the above eqn becomes, 

i N 3 
[_y Vi 2 -4 y[ e -«(r«-ro)( e -a(r«-ro) _ 2 )1 

L 2^ se ar °(e ar ° - 16) ^ L v ;J 

2 i=i 

= -/io^(r) (39) 
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1 N 

[L ]T Vi 2 - 7 ^[ e - a ^- r °)(e- a( ^- ro) - 2)] 

i=l i<j 
1 ^ 

--E[^ 2 + ^ 2 + A ^ 2 M f ) 

= -/io^(r) (40) 

Now for a = .29 and r = 9.758 (both in oscillator units) [32]. We have 
checked that for these choice of parameters, Morse solution is extremely 
good, a = 52 x 10~ 10 cm, s = 12 x 10~ 7 cm, the interaction strength 7 is 
given by 

= 4 o« = 2M x 1Q -5 / 41 x 

se Qr °(e ar ° - 16) 

For mean field calculation the value of interaction strength was taken to 
be 4.33 x 10~ 3 . For this problem we are interested in the limit 7 << 1. 
The case 7 >> 1 is usually known as the Thomas Fermi limit. For 7 = 
2.64 x 10~ 5 , the eigenvalue equation reduces to a minimally perturbed system 
of d dimensional anisotropic oscillator where d = 3N and N is the number of 
particles. The whole concept of bound states of Morse dimers is outside the 
range of this limit, so the nonexistence of two-body bound states is ensured 
by choosing the above parameters. 

Even though 7 << 1, we solve the eigenvalue eqn nonperturba- 
tively with Generalized Feynman-Kac procedure. Energies and frequencies at 
zero temperature are obtained by solving Eq. (4) and using Eq.(ll). To cal- 
culate the analogous quantities at finite temperature we use Eq.(20). We can 
get the energy of both condensate and noncondensate using Eq.(4) Eq.(20 ) 
by generating a large number of paths and then averaging the results for all 
the paths. Since original Feynman-Kac method [14,15] is computationally 
inefficient we incorporate importance sampling in our random walk and use 
trial function of the form given in Eq.(26) 

Evaluation of temperature dependent mode frequencies 
: Following the prescription in [4], we see that for a fixed 'N' relationship 
H(N ,T) ~ fi(N,T = 0)(f ) 2 / 5 or equivalent^ N /N = [^§Sj] 5/2 gener- 
ates the condensation fraction N /N as a function of time. One can generate 
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this from experiments also. From the thermodynamical limit we get iVo as a 
function of time and run our zero temperature code with same number of N 
as the dynamics of the finite temperature condensate are essentially the same 
as those of a zero temperature condensate with the same value of N . This 
is called method I. The other way to calculate energy is to run the code for 
different simulation times corresponding to different temperatures as time is 
defined as inverse temperature. This is identified as method II. /i(Nq, T) and 
fi{N ,T = 0) are calculated using fi = jj{fikin + Hho + 2fx int ) and fi k in, Hho 
and Hi n t are calculated as described in Ref[25]. Later in Fig. 5 (our data) of 
Section 4.2, we see the effects of interaction on the condensation fraction. 
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4 Results 



4.1 Excitation spectra at T=0 for different symme- 
tries 

We chose Rb 87 as an example of a weakly interacting dilute Bose gas as in 
Ref[29]. We simulate 2000 Rb atoms interacting via the Morse potential. We 
choose a = 52 x 10 _10 cm, length scale s of the problem as 12 x 10~ 7 cm. In the 
following table, we show the ground state energy of the particle. With repul- 
sive interactions, the Energy/particle increases with the incresae in number 
of particles in the trap [Fig. 1-2] whereas the energy gap between the different 
symmetry states decreases as evident from Fig. 3. However we see a differ- 
ent trend [Fig. 4] for m = mode where the excitation frequencies increase 
with the increase in number of atoms. This agrees with the Hartree-Fock 
spectrum in the Fig 2 of Ref[33]. 

Table 1: Results for the ground state energy of 2000 Rb 87 atoms in 
a trap with uj x — u y — 1, A = uj z — v^8 in the interacting case ; The 
table shows how energy varies with the number of particles in the 
Gross Pitaevski(GP) case [34] and GFK method. 



N 


E/N(GP) 


E/N(GFK) 


1 


2.414 


2.414 213 


100 


2.66 


2.52230(5) 


200 


2.86 


2.63141(1) 


500 


3.30 


2.9588(4) 


1000 


3.84 


3.5047(4) 


2000 


4.61 


4.5962(3) 
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Number of Atoms [N] 



Figure 1: A plot for the Condensate Energy/Atom versus Number of atoms 
in trap for 2000 oarticles for the ground state ; this work 




Number of Atoms [N] 

Figure 2: A plot for the Condensate Energy/Particle versus Number of atoms 
in trap for 2000 particles for the m = 2 mode; this work 
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Table 2: frequency u for lowest lying modes 



N 


mode of oscillation 


energy 


a; (this work) 


a;(JILA TOP) 


2000 


ground state 


4.596(3) 






2000 


m = 2 


5.860(1) 


1.264(4) 


1.4 


2000 


m = 


7.295(3) 


2.699(6) 


1.8 




Number of Atoms [N] 

Figure 3: A plot of Excitation Frequency vs Number of Atoms for lowest 
lying m = 2 mode ;this work 
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500 1000 1500 

Number of Atoms [N] 



Figure 4: A plot of Excitation Frequency vs Number of Atoms for lowest 
lying m = mode ; this work 

4.2 Effects of temperature on condensation fraction 

Density of condensate atoms decreases in the trap as temperature increases. 
This lowers the interaction energy of the condensate atoms resulting in a 
shift in the critical temperature. As a matter of fact in the interacting case, 
the critical temperature decreases. This is a very unique feature of trapped 
gas. In the case of uniform gas we see an opposite trend. 
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Reduced Temperature [ T/T () ] 

Figure 5: Condensation fraction vs Reduced Temperature ; this work. The 
middle curve corresponds to the 2000 interacting atoms(GFK sumlation 
data;ourwork) and the outer one corresponds to the noninteracting case. 
The innermost curve also corresponds interacting atoms[Ref 37]. The number 
of condensed particles decreases with the interaction 

4.3 Effects of temperature on the frequency shifts; 
comparison with other experiments and theories 

Next we give an account of how our path integral simulations compare with 
other theoretical and experimental data. Fig 3a in Ref[2] represents JILA 
TOP data where one observes a large temperature dependent frequency shift 
for both m=0 and m=2 modes. For m=2 mode, starting from Stringari 
limit it decreases all the way up to 0.9T whereas for m=0 mode it shows 
a rising trend with rise in temperature. Our path integral data for m = 2 
mode in Fig [6] shows similar decreasing trend as JILA TOP data in Fig 3a 
of Ref[2] and best theoretical data in Fig 1 of Ref[7] all the way up to 0.9T 
whereas data generated by Hartree-Fock-Bogoliubov[HFB] method in Fig 1 
of Ref[4] agree with JILA TOP only up to 0.7To . Finally our data Fig[7] 
for temperature variation of m=0 mode agrees with JILA data in Fig[3a] of 
Ref[2] and Morgan data in Fig 1 of Ref[7] when the effect of thermal cloud 
is considered. 
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Figure 6: Effects of temperature on m = 2 mode; this work 
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Reduced Temperature [ T/T Q ] 

Figure 7: Effects of temperature onm = mode from GFK considering non- 
condensate dynamicsfthis work], shows resemblence with JILA[2] and Mor- 
gan et al[7] 

4.4 Discussions 

We are solving the full Hamiltonian with some realistic potential and thereby 
trying to solve the many body problem fully quantum mechanically and non- 
perturbatively. So for both m = and m = 2 modes, our calculations adopt- 
ing Feynman Kac path integral technique represent the collective behavior 
of the Bose gas. Our work [Fig. 6] agrees with JILA experiment [2] for m=2 
mode. The other theoretical work shows the reverse trend [Fig 1 of Ref 4]. 
We found that considering the dynamics of condensates alone and the effect 
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of finite temperature as static thermal cloud, we do not achieve the upward 
shifts of frequencies as shown by JILA data for m=0 mode. In fact, we agree 
with Ref [5] that as m=0 happens to be a coupled mode, we need to consider 
the dynamics of thermal cloud to obtain a satisfactory agreemenent with the 
experimental data. Eventually We consider the dynamics of thermal cloud 
and get the upward shift of data [Fig 7] as observed in JILA [2] and Ref [7, 8]. 

At T=0 for m = 2, we observe that as N increases the energies 
grow, but the splitting between the ground and excited state decreases - an 
essential feature of Bose Condensation. But in the similar case for m = 
mode we observe that both the energies and the gap between the ground and 
higher excited states increase with increase in number of atoms. This agrees 
with Hartree-Fock spectrum of Ref [33] where the author had to improve these 
results by using random phase approximation to get an agreement with the 
experimental results. For m = 2 mode the Stringari limit turns out to be 
1.264(4) as opposed to the experimental value of 1.4. On the other hand 
m = mode the Stringari limit is 2.699(6) as opposed to the experimental 
value of 1.8. So for the coupled m = mode, our path integral method 
does not work better than Hartree-Fock theory and thereby does not yield a 
correct value of Stringari limit. The reason that the quantitative agreement 
between the Stringari limit predicted by us particularly for m = mode and 
the experiments is not good, might be the use of Gaussian wave functions 
as the trial functions. In general, the frequencies of collective modes do not 
have direct correspondence to harmonic oscillator states because they do not 
include correlations. It is legitimate to use harmonic oscillator solutions as 
trial functions for m=2 mode as JILA experiment does not correspond to the 
Thomas Fermi limit [35]. But for m=0 mode we need to include correlations 
in the wave function as m=0 is a coupled mode. nthebibliography99 
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5 Conclusions: 



We have used GFK to bring out the many body effects between the cold 
Rb atoms. Numerical work with bare Feynman-Kac procedure employing 
modern computers was reported[15] for the first time for few electron systems 
after forty years of the original work[14] and seemed to be really useful for 
calculating atomic ground states[19]. A fairly good success in atomic physics 
motivated us to apply it to Condensed matter Physics. 

Gross-Pitaevskii (GP) technique [33] does not include correlations 
in the solutions explicitly and calculates energy at the variational level. These 
energies are upper bounds to the actual energies of the system. We have been 
successful in achieving a lower value for Rb ground state than that obtained 
by GP. This correct trend in our calculated energies for different symmetry 
states enables one to calculate the frequencies more accurately by Feynman 
Kac path integral method. As a result, Diffusion Monte Carlo codes based 
on nonperturbative quantum appraoach can handle temperature very accu- 
rately and we do not see any breakdown near T c . For the first time we have 
calculated finite temperature properties beyond mean field approximation 
by Quantum Monte Carlo technique. The only other non mean field calcula- 
tions at T = worth mentioning in this context is the work done by Blume 
et al[36]. We have calculated spectrum of Rb gas by considering realistic 
potentials like Morse potential etc. instead of conventional pseudopotentials 
for the first time. 

We have been able to calculate the lowest lying excitation frequen- 
cies for m = and m = 2 modes by Feynman-Kac path integral technique 
in a very simple way. We have found an alternative to Gross-Pitaevskii tech- 
nique and other mean field calculations which works at all the temperature. 
Simulating 2000 atoms with the path intgral method we have been able to 
capture some of the signatures of Bose condensation like decrease of excita- 
tion frequencies with number of atoms, lowering of condensation fraction in 
the interacting case etc. In our non mean field study, we see agreement with 
experimental study all the way to T = 0.9T C [Fig. 6]. This is because of 
the fact that we have been able to solve the related many body theory very 
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accurately with the nonperturbative and quantum mechanical approach. At 
this point our results agree with the experimental results only qualitatively 
as we are restricting ourselves to the choice of Gaussin trial functions. To 
improve our results quantitatively we need to use correlated trial functions. 
This would be a nontrivial extention of the present work and will be reported 
elsewhere. The simplicity in our method is appealing as it is extremely easy 
to implement and our fortran code at this point consists of about 270 lines. 
In fact mere ability to add, subtract and toss a coin enables one to solve 
many body theory with our path integral technique. 

We employ an algorithm which is essentially parallel in nature 
so that eventually we can parallelize our code and calculate thermodynamic 
properties of bigger systems taking advantage of new computer architech- 
tures. This work is in progress. We are continuing on this problem and hope 
that this technique will inspire others to do similar calculations. 
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